#bbox_of_interest = [-59.958, -9.818, -58.633 , -8.946 ]
#search = catalog.search(collections=["io-biodiversity"], bbox=bbox_of_interest)
#item = next(search.get_items())
#item
everest = [86.9250, 27.9881]
seattle = [-122.332, 47.6062]
grand_canyon = [-112.107676, 36.101690]
mount_fuji = [138.7274, 35.3606]
mont_blanc = [6.865000, 45.832778]
imjalake = [ 86.92586485793316, 27.899354509705653]
areas_of_interestAlaska = {
"type": "Polygon",
"coordinates": [
[
[-148.56536865234375, 60.80072385643073],
[-147.44338989257812, 60.80072385643073],
[-147.44338989257812, 61.18363894915102],
[-148.56536865234375, 61.18363894915102],
[-148.56536865234375, 60.80072385643073],
]
],
}
areas_of_interestImjalake = {
"type": "Polygon",
"coordinates": [
[
[
86.92228070099549,
27.9096329170643
],
[
86.92228070099549,
27.882026779073712
],
[
86.95880774284717,
27.882026779073712
],
[
86.95880774284717,
27.9096329170643
],
[
86.92228070099549,
27.9096329170643
]
]
],
}
areas_of_interestLhonaklake = {
"type": "Polygon",
"coordinates": [
[
[
88.31558176501454,
27.97312474305086
],
[
88.31558176501454,
27.937831970683348
],
[
88.34995817854201,
27.937831970683348
],
[
88.34995817854201,
27.97312474305086
],
[
88.31558176501454,
27.97312474305086
]
]
],
}
areas_of_interestLhonaklaketoDam = {
"type": "Polygon",
"coordinates": [
[
[
88.31217526255193,
27.960400536405857
],
[
88.31217526255193,
27.744973669773245
],
[
88.57021256641684,
27.744973669773245
],
[
88.57021256641684,
27.960400536405857
],
[
88.31217526255193,
27.960400536405857
]
]
],
}
areas_of_interestImjalakeold = {
"type": "Polygon",
"coordinates": [
[
[
86.90840990978438,
27.914442874882738
],
[
86.90840990978438,
27.87583921087139
],
[
86.96126773837943,
27.87583921087139
],
[
86.96126773837943,
27.914442874882738
],
[
86.90840990978438,
27.914442874882738
]
]
],
}
areas_of_interestImjalakeGL = {
"type": "Polygon",
"coordinates": [
[
[
88.14439581658235,
27.929979855917935
],
[
88.14439581658235,
27.90157666024274
],
[
88.21527228069135,
27.90157666024274
],
[
88.21527228069135,
27.929979855917935
],
[
88.14439581658235,
27.929979855917935
]
]
],
}
#areas_of_interest = {"type": "Point", "coordinates": imjalake}
#areas_of_interest = areas_of_interestImjalake
#areas_of_interest = areas_of_interestLhonaklake
areas_of_interest = areas_of_interestImjalake
time_of_interest = "2019-06-01/2023-10-01"
time_of_interest2 = "2021-08-01/2022-12-01"
yearofinterest = str(2021)
time_of_interest = yearofinterest + "-07-01/" + yearofinterest + "-08-01"
time_of_interest
'2021-07-01/2021-08-01'
import rasterio
from rasterio import windows
from rasterio import features
from rasterio import warp
import numpy as np
from PIL import Image
def displayRawImage(least_cloudy_item):
asset_href = least_cloudy_item.assets["visual"].href
with rasterio.open(asset_href) as ds:
aoi_bounds = features.bounds(areas_of_interest)
print("area of interest")
print(aoi_bounds)
warped_aoi_bounds = warp.transform_bounds("epsg:4326", ds.crs, *aoi_bounds)
print("wraped aoi bounds")
# print(warped_aoi_bounds)
aoi_window = windows.from_bounds(transform=ds.transform, *warped_aoi_bounds)
print(aoi_window)
band_data = ds.read(window=aoi_window)
#print(band_data)
img = Image.fromarray(np.transpose(band_data, axes=[1, 2, 0]))
w = img.size[0]
h = img.size[1]
aspect = w / h
target_w = 800
target_h = (int)(target_w / aspect)
return img.resize((target_w, target_h), Image.Resampling.BILINEAR)
print(time_of_interest2)
2021-08-01/2022-12-01
def displayRawImagesentinel1(least_cloudy_item):
img_arr = rioxarray.open_rasterio(least_cloudy_item.assets["rendered_preview"].href).squeeze()
#print(item.datetime.strftime("%x"))
img_arr.plot.imshow(figsize=(25, 20));
return
from IPython.display import Image
from pystac.extensions.eo import EOExtension as eo
import pystac_client
import planetary_computer
import pystac_client
import planetary_computer
catalog = pystac_client.Client.open(
"https://planetarycomputer.microsoft.com/api/stac/v1/",
modifier=planetary_computer.sign_inplace,
)
import numpy as np
from PIL import Image
import rioxarray
import matplotlib.pyplot as plt
def displayRawImage222222(least_cloudy_item):
asset_href = least_cloudy_item.assets["visual"].href
with rasterio.open(asset_href) as ds:
aoi_bounds = features.bounds(areas_of_interest)
print("area of interest")
print(aoi_bounds)
warped_aoi_bounds = warp.transform_bounds("epsg:4326", ds.crs, *aoi_bounds)
print("wraped aoi bounds")
# print(warped_aoi_bounds)
aoi_window = windows.from_bounds(transform=ds.transform, *warped_aoi_bounds)
print(aoi_window)
band_data = ds.read(window=aoi_window)
#print(band_data)
return band_data[2]
def displayRawImageSAC(least_cloudy_item):
asset_href = least_cloudy_item.assets["SCL"].href
with rasterio.open(asset_href) as ds:
aoi_bounds = features.bounds(areas_of_interest)
print("area of interest")
print(aoi_bounds)
warped_aoi_bounds = warp.transform_bounds("epsg:4326", ds.crs, *aoi_bounds)
print("wraped aoi bounds")
# print(warped_aoi_bounds)
aoi_window = windows.from_bounds(transform=ds.transform, *warped_aoi_bounds)
print(aoi_window)
band_data = ds.read(window=aoi_window)
#print(band_data)
band_datashow = band_data == 6
# plt.imshow(band_datashow[0])
return band_data[0]
for i in range(2012, 2024, 1):
print(i)
yearofinterest = str(i)
time_of_interest = yearofinterest + "-08-01/" + yearofinterest + "-11-01"
print(time_of_interest)
search = catalog.search(collections=["sentinel-2-l2a"], intersects=areas_of_interest, datetime=time_of_interest, query={"eo:cloud_cover": {"lt": 50, "gt":2 }})
# search = catalog.search(collections=["sentinel-2-l2a"], intersects=areas_of_interest, datetime=time_of_interest )
items = search.item_collection()
print(f"Returned {len(items)} Items")
if (len(items) >0):
least_cloudy_item = min(items, key=lambda item: eo.ext(item).cloud_cover)
item = least_cloudy_item
answer = displayRawImage(least_cloudy_item)
plt.imshow(answer)
plt.title(item.datetime.strftime("%B %d, %Y"))
plt.show()
answer = displayRawImageSAC(least_cloudy_item)
print(item.datetime.strftime("%x"))
plt.imshow(answer)
plt.title(item.datetime.strftime("%B %d, %Y"))
plt.show()
2012 2012-08-01/2012-11-01 Returned 0 Items 2013 2013-08-01/2013-11-01 Returned 0 Items 2014 2014-08-01/2014-11-01 Returned 0 Items 2015 2015-08-01/2015-11-01 Returned 0 Items 2016 2016-08-01/2016-11-01 Returned 2 Items area of interest (86.92228070099549, 27.882026779073712, 86.95880774284717, 27.9096329170643) wraped aoi bounds Window(col_off=9239.013502126181, row_off=1282.556031809363, width=359.63644153829955, height=305.9773380664992)
area of interest (86.92228070099549, 27.882026779073712, 86.95880774284717, 27.9096329170643) wraped aoi bounds Window(col_off=4619.506751063091, row_off=641.2780159046815, width=179.81822076914978, height=152.9886690332496) 10/30/16
2017 2017-08-01/2017-11-01 Returned 3 Items area of interest (86.92228070099549, 27.882026779073712, 86.95880774284717, 27.9096329170643) wraped aoi bounds Window(col_off=9239.013502126181, row_off=1282.556031809363, width=359.63644153829955, height=305.9773380664992)
area of interest (86.92228070099549, 27.882026779073712, 86.95880774284717, 27.9096329170643) wraped aoi bounds Window(col_off=4619.506751063091, row_off=641.2780159046815, width=179.81822076914978, height=152.9886690332496) 10/30/17
2018 2018-08-01/2018-11-01 Returned 8 Items area of interest (86.92228070099549, 27.882026779073712, 86.95880774284717, 27.9096329170643) wraped aoi bounds Window(col_off=9239.013502126181, row_off=1282.556031809363, width=359.63644153829955, height=305.9773380664992)
area of interest (86.92228070099549, 27.882026779073712, 86.95880774284717, 27.9096329170643) wraped aoi bounds Window(col_off=4619.506751063091, row_off=641.2780159046815, width=179.81822076914978, height=152.9886690332496) 10/30/18
2019 2019-08-01/2019-11-01 Returned 6 Items area of interest (86.92228070099549, 27.882026779073712, 86.95880774284717, 27.9096329170643) wraped aoi bounds Window(col_off=9239.013502126181, row_off=1282.556031809363, width=359.63644153829955, height=305.9773380664992)
area of interest (86.92228070099549, 27.882026779073712, 86.95880774284717, 27.9096329170643) wraped aoi bounds Window(col_off=4619.506751063091, row_off=641.2780159046815, width=179.81822076914978, height=152.9886690332496) 10/15/19
2020 2020-08-01/2020-11-01 Returned 7 Items area of interest (86.92228070099549, 27.882026779073712, 86.95880774284717, 27.9096329170643) wraped aoi bounds Window(col_off=9239.013502126181, row_off=1282.556031809363, width=359.63644153829955, height=305.9773380664992)
area of interest (86.92228070099549, 27.882026779073712, 86.95880774284717, 27.9096329170643) wraped aoi bounds Window(col_off=4619.506751063091, row_off=641.2780159046815, width=179.81822076914978, height=152.9886690332496) 10/29/20
2021 2021-08-01/2021-11-01 Returned 2 Items area of interest (86.92228070099549, 27.882026779073712, 86.95880774284717, 27.9096329170643) wraped aoi bounds Window(col_off=9239.013502126181, row_off=1282.556031809363, width=359.63644153829955, height=305.9773380664992)
area of interest (86.92228070099549, 27.882026779073712, 86.95880774284717, 27.9096329170643) wraped aoi bounds Window(col_off=4619.506751063091, row_off=641.2780159046815, width=179.81822076914978, height=152.9886690332496) 10/14/21
2022 2022-08-01/2022-11-01 Returned 5 Items area of interest (86.92228070099549, 27.882026779073712, 86.95880774284717, 27.9096329170643) wraped aoi bounds Window(col_off=9239.013502126181, row_off=1282.556031809363, width=359.63644153829955, height=305.9773380664992)
area of interest (86.92228070099549, 27.882026779073712, 86.95880774284717, 27.9096329170643) wraped aoi bounds Window(col_off=4619.506751063091, row_off=641.2780159046815, width=179.81822076914978, height=152.9886690332496) 10/24/22
2023 2023-08-01/2023-11-01 Returned 2 Items area of interest (86.92228070099549, 27.882026779073712, 86.95880774284717, 27.9096329170643) wraped aoi bounds Window(col_off=9239.013502126181, row_off=1282.556031809363, width=359.63644153829955, height=305.9773380664992)
area of interest (86.92228070099549, 27.882026779073712, 86.95880774284717, 27.9096329170643) wraped aoi bounds Window(col_off=4619.506751063091, row_off=641.2780159046815, width=179.81822076914978, height=152.9886690332496) 09/29/23
def displayRawImageSNOW(least_cloudy_item):
img_arr = rioxarray.open_rasterio(least_cloudy_item.assets["visual"].href).squeeze()
print(item.datetime.strftime("%x"))
img_arr.plot.imshow(figsize=(25, 20));
return
def displayRawImageSNOW2(least_cloudy_item):
img_arr = rioxarray.open_rasterio(least_cloudy_item.assets["visual"].href).squeeze()
print(item.datetime.strftime("%x"))
return img_arr
def displayRawImageSNOWCLASSIFY(least_cloudy_item):
img_arr = rioxarray.open_rasterio(least_cloudy_item.assets["SCL"].href).squeeze()
# print(least_cloudy_item.datetime.strftime("%x"))
#https://custom-scripts.sentinel-hub.com/custom-scripts/sentinel-2/scene-classification/
img_arrnew = img_arr == 11 #11-snow, 6-water
return img_arrnew
def displayRawImageSNOWCLASSIFY2(least_cloudy_item):
img_arr = rioxarray.open_rasterio(least_cloudy_item.assets["SCL"].href).squeeze()
# print(least_cloudy_item.datetime.strftime("%x"))
#https://custom-scripts.sentinel-hub.com/custom-scripts/sentinel-2/scene-classification/
img_arrnew = img_arr == 11 #11-snow, 6-water
return img_arr
from datashader.transfer_functions import shade, stack
from datashader.colors import Elevation
from xrspatial import hillshade
import pystac_client
import planetary_computer
import rioxarray
search = catalog.search(collections=["nasadem"], intersects=areas_of_interest)
items = search.item_collection()
print(f"Returned {len(items)} Items")
item = items[0]
da = (
rioxarray.open_rasterio(item.assets["elevation"].href)
.squeeze()
.drop("band")[:-1, :-1]
.coarsen({"y": 5, "x": 5})
.mean()
)
Returned 1 Items
from datashader.transfer_functions import shade
import matplotlib.pyplot as plt
from xrspatial.classify import natural_breaks
for i in range(1):
print(i)
yearofinterest = str(i)
time_of_interest = yearofinterest + "-06-01/" + yearofinterest + "-08-01"
print(time_of_interest)
#search = catalog.search(collections=["nasadem"], intersects=areas_of_interest, query={"eo:cloud_cover": {"lt": 450, "gt":2 }})
search = catalog.search(collections=["nasadem"], intersects=areas_of_interest )
items = search.item_collection()
print(f"Returned {len(items)} Items")
if (len(items) >0):
item = items[0]
answer = rioxarray.open_rasterio(item.assets["elevation"].href).squeeze()
natural_breaks_result = natural_breaks(answer, num_sample=20000, k=15)
shade(natural_breaks_result, cmap=plt.get_cmap("terrain"), how="linear")
# classify the image using natural breaks
0 0-06-01/0-08-01 Returned 1 Items
least_cloudy_item = min(items, key=lambda item: eo.ext(item).cloud_cover)
print(
f"Choosing {least_cloudy_item.id} from {least_cloudy_item.datetime.date()}"
f" with {eo.ext(least_cloudy_item).cloud_cover}% cloud cover"
)
item = least_cloudy_item
asset_href = least_cloudy_item.assets["visual"].href
items = list(search.items())
for item in items:
print(item)
import rich.table
asset_table = rich.table.Table("Asset Key", "Asset Title")
for key, value in items[-1].assets.items():
asset_table.add_row(key, value.title)
asset_table
import rasterio
from rasterio import windows
from rasterio import features
from rasterio import warp
import numpy as np
from PIL import Image
asset_href = least_cloudy_item.assets["visual"].href
with rasterio.open(asset_href) as ds:
aoi_bounds = features.bounds(areas_of_interest)
print("area of interest")
print(aoi_bounds)
warped_aoi_bounds = warp.transform_bounds("epsg:4326", ds.crs, *aoi_bounds)
print("wraped aoi bounds")
print(warped_aoi_bounds)
aoi_window = windows.from_bounds(transform=ds.transform, *warped_aoi_bounds)
print(aoi_window)
band_data = ds.read(window=aoi_window)
print(band_data)
area of interest (88.14439581658235, 27.90157666024274, 88.21527228069135, 27.929979855917935) wraped aoi bounds (612596.7648950818, 3086825.9764255327, 619602.1054647851, 3090039.936055328) Window(col_off=1259.6764895081797, row_off=998.006394467142, width=700.534056970333, height=321.39596297958633) [[[198 255 255 ... 182 193 200] [212 255 255 ... 189 202 202] [245 255 255 ... 197 208 206] ... [255 255 255 ... 64 69 59] [255 255 255 ... 75 72 60] [255 255 255 ... 80 75 69]] [[204 255 255 ... 161 168 171] [182 255 255 ... 166 177 176] [171 255 255 ... 174 184 181] ... [255 255 255 ... 69 67 58] [255 255 255 ... 83 69 64] [255 255 255 ... 82 80 94]] [[166 255 255 ... 117 124 125] [141 255 255 ... 120 129 127] [113 255 255 ... 126 135 132] ... [255 255 255 ... 73 65 58] [255 255 255 ... 77 78 87] [255 255 255 ... 74 86 94]]]
areas_of_interest
{'type': 'Polygon',
'coordinates': [[[88.14439581658235, 27.929979855917935],
[88.14439581658235, 27.90157666024274],
[88.21527228069135, 27.90157666024274],
[88.21527228069135, 27.929979855917935],
[88.14439581658235, 27.929979855917935]]]}
import rasterio
from rasterio import windows
from rasterio import features
from rasterio import warp
import numpy as np
from PIL import Image
asset_href = least_cloudy_item.assets["SCL"].href
with rasterio.open(asset_href) as ds:
aoi_bounds = features.bounds(areas_of_interest)
print(aoi_bounds)
warped_aoi_bounds = warp.transform_bounds("epsg:4326", ds.crs, *aoi_bounds)
print(warped_aoi_bounds)
aoi_window = windows.from_bounds(transform=ds.transform, *warped_aoi_bounds)
print(aoi_window)
band_data = ds.read(window=aoi_window)
print(band_data)
(88.14439581658235, 27.90157666024274, 88.21527228069135, 27.929979855917935) (612596.7648950818, 3086825.9764255327, 619602.1054647851, 3090039.936055328) Window(col_off=629.8382447540898, row_off=499.003197233571, width=350.2670284851665, height=160.69798148979316) [[[11 11 11 ... 5 5 5] [11 11 11 ... 5 5 5] [11 11 11 ... 5 5 5] ... [11 11 11 ... 7 2 2] [11 11 11 ... 7 7 2] [11 11 11 ... 7 7 7]]]
band_data.shape
(1, 161, 350)
type(band_data)
numpy.ndarray
This item's data can be read in using xarray's open_rasterio function to load the data into an array:
import rioxarray
item = least_cloudy_item
img_arr = rioxarray.open_rasterio(item.assets["visual"].href).squeeze()
print(item.datetime.strftime("%x"))
img_arr.plot.imshow(figsize=(25, 20));
08/12/22
img_arr
<xarray.DataArray (y: 5490, x: 5490)>
array([[11, 11, 11, ..., 5, 5, 5],
[11, 11, 11, ..., 5, 5, 5],
[11, 11, 11, ..., 5, 5, 5],
...,
[ 9, 9, 9, ..., 4, 4, 4],
[ 9, 9, 9, ..., 4, 4, 4],
[ 9, 9, 9, ..., 4, 4, 4]], dtype=uint8)
Coordinates:
band int64 1
* x (x) float64 6e+05 6e+05 6e+05 ... 7.098e+05 7.098e+05 7.098e+05
* y (y) float64 3.1e+06 3.1e+06 3.1e+06 ... 2.99e+06 2.99e+06
spatial_ref int64 0
Attributes:
AREA_OR_POINT: Area
_FillValue: 0
scale_factor: 1.0
add_offset: 0.0
img_arrnew = img_arr == 11 #11-snow, 6-water
img_arrnew
<xarray.DataArray (y: 5490, x: 5490)>
array([[ True, True, True, ..., False, False, False],
[ True, True, True, ..., False, False, False],
[ True, True, True, ..., False, False, False],
...,
[False, False, False, ..., False, False, False],
[False, False, False, ..., False, False, False],
[False, False, False, ..., False, False, False]])
Coordinates:
band int64 1
* x (x) float64 6e+05 6e+05 6e+05 ... 7.098e+05 7.098e+05 7.098e+05
* y (y) float64 3.1e+06 3.1e+06 3.1e+06 ... 2.99e+06 2.99e+06
spatial_ref int64 0snowcount = sum((sum(img_arrnew)))
print("Snow area = " + str(snowcount))
Snow area = <xarray.DataArray ()>
array(2291118)
Coordinates:
band int64 1
spatial_ref int64 0
type(snowcount)
xarray.core.dataarray.DataArray
snowcount.values
array(2291118)
type(img_arrnew)
xarray.core.dataarray.DataArray
import rioxarray
img_arr = rioxarray.open_rasterio(item.assets["elevation"].href).squeeze().drop("band")
img_arr.plot.imshow(figsize=(15, 10));
Y